###################################################
### chunk number 1: setup
###################################################
#line 47 "predictionet_pqg_short_course"
library(pgfSweave)
setCacheDir("cache")
options(keep.source=TRUE)


###################################################
### chunk number 2: initseed
###################################################
#line 53 "predictionet_pqg_short_course"
set.seed(123)


###################################################
### chunk number 3: preparedata
###################################################
#line 130 "predictionet_pqg_short_course"
## read the gene list with the corresponding annotations
pik3ca.list <- read.csv("files/loi2010_pik3cags_signature_278.csv", stringsAsFactors=FALSE)
rownames(pik3ca.list) <- pik3ca.list[ , "probe"]
## remove special characters in the gene symbols
nn <- gsub(pattern="[-]|[+]|[*]|[%]|[$]|[#]|[{]|[}]|[[]|[]]|[|]|[\\^]", replacement=".", x=pik3ca.list[ ,"Gene.symbol"])
pik3ca.list <- cbind(pik3ca.list, "Gene.symbol.curated"=toupper(nn))

## remove unknown and duplicated genes by respecting the order of the list
pik3ca.list <- pik3ca.list[is.na(pik3ca.list[ ,"Gene.symbol.curated"]) | !duplicated(pik3ca.list[ ,"Gene.symbol.curated"]), , drop=FALSE]
## select only the top genes
genen <- 20
if(genen < 2 || genen > nrow(pik3ca.list)) { genen <- nrow(pik3ca.list) }
pik3ca.list <- pik3ca.list[1:genen, , drop=FALSE]

## load the first dataset
load("files/marray_breast1/minn2007_frma.RData")
data1 <- data[ ,rownames(pik3ca.list), drop=FALSE]
annot1 <- annot[rownames(pik3ca.list), , drop=FALSE]
colnames(data1) <- rownames(annot1) <- pik3ca.list[ ,"Gene.symbol.curated"]
demo1 <- demo
## save the data in a CSV file
write.csv(t(data1), file="files/marray_breast1_pik3ca.csv")
## free memory
rm(list=c("data", "data.bin", "annot", "demo"))
gc()

## load the second dataset
load("files/marray_breast2/transbig2006affy_frma.RData")
data2 <- data[ ,rownames(pik3ca.list), drop=FALSE]
annot2 <- annot[rownames(pik3ca.list), , drop=FALSE]
colnames(data2) <- rownames(annot2) <- pik3ca.list[ ,"Gene.symbol.curated"]
demo2 <- demo
## save the data in a CSV file
write.csv(t(data2), file="files/marray_breast2_pik3ca.csv")
## free memory
rm(list=c("data", "data.bin", "annot", "demo"))
gc()


###################################################
### chunk number 4: preparepriors
###################################################
#line 182 "predictionet_pqg_short_course"
## read priors generated by the Predictive Networks web application
pn.priors <- read.csv("files/priors_pik3ca.csv", stringsAsFactors=FALSE)
## the column names should be: subject, predicate, object, direction, evidence, sentence, article, network

## remove special characters in the gene symbols
pn.priors[ ,"subject"] <- gsub(pattern="[-]|[+]|[*]|[%]|[$]|[#]|[{]|[}]|[[]|[]]|[|]|[\\^]", replacement=".", x=pn.priors[ ,"subject"])
pn.priors[ ,"object"] <- gsub(pattern="[-]|[+]|[*]|[%]|[$]|[#]|[{]|[}]|[[]|[]]|[|]|[\\^]", replacement=".", x=pn.priors[ ,"object"])

## missing values
pn.priors[!is.na(pn.priors) & (pn.priors == "" | pn.priors == " " | pn.priors == "N/A")] <- NA

## select only the interactions in which the genes are comprised in our gene expression dataset
myx <- is.element(pn.priors[ , "subject"], pik3ca.list[ ,"Gene.symbol.curated"]) & is.element(pn.priors[ , "object"], pik3ca.list[ ,"Gene.symbol.curated"])
pn.priors <- pn.priors[myx, , drop=FALSE]


###################################################
### chunk number 5: tablepriors
###################################################
#line 201 "predictionet_pqg_short_course"
## print the first rows of 'pn.priors'
print(head(pn.priors))


###################################################
### chunk number 6: preparepriorscount
###################################################
#line 208 "predictionet_pqg_short_course"
## build prior counts
priors.count <- matrix(0, nrow=nrow(pik3ca.list), ncol=nrow(pik3ca.list), dimnames=list(pik3ca.list[ ,"Gene.symbol.curated"], pik3ca.list[ ,"Gene.symbol.curated"]))

for(i in 1:nrow(pn.priors)) {
	switch(tolower(pn.priors[i, "direction"]), 
	"right"={ priors.count[pn.priors[i, "subject"], pn.priors[i, "object"]] <- priors.count[pn.priors[i, "subject"], pn.priors[i, "object"]] + ifelse(!is.na(pn.priors[i, "evidence"]), ifelse(tolower(pn.priors[i, "evidence"]) == "positive", +1, -1), 0) }, 
	"left"={ priors.count[pn.priors[i, "object"], pn.priors[i, "subject"]] <- priors.count[pn.priors[i, "object"], pn.priors[i, "subject"]] +  ifelse(!is.na(pn.priors[i, "evidence"]), ifelse(tolower(pn.priors[i, "evidence"]) == "positive", +1, -1), 0) }, 
	{ priors.count[pn.priors[i, "subject"], pn.priors[i, "object"]] <- priors.count[pn.priors[i, "subject"], pn.priors[i, "object"]] +  ifelse(!is.na(pn.priors[i, "evidence"]), ifelse(tolower(pn.priors[i, "evidence"]) == "positive", +1, -1), 0)
	if(pn.priors[i, "object"] != pn.priors[i, "subject"]) { priors.count[pn.priors[i, "object"], pn.priors[i, "subject"]] <- priors.count[pn.priors[i, "object"], pn.priors[i, "subject"]] +  ifelse(!is.na(pn.priors[i, "evidence"]), ifelse(tolower(pn.priors[i, "evidence"]) == "positive", +1, -1), 0) } })
}
## negative count represent evidence for ABSENCE of an interaction, positive otherwise


###################################################
### chunk number 7: priorscounttab
###################################################
#line 224 "predictionet_pqg_short_course"
print(table(priors.count))


###################################################
### chunk number 8: installpackage eval=FALSE
###################################################
## #line 314 "predictionet_pqg_short_course"
## install.packages(pkgs="predictionet_0.7.tar.gz", repos=NULL)


###################################################
### chunk number 9: installpackagedep eval=FALSE
###################################################
## #line 320 "predictionet_pqg_short_course"
## install.packages("penalized", repos="http://cran.r-project.org", dependencies=TRUE)


###################################################
### chunk number 10: loadpredictionet
###################################################
#line 329 "predictionet_pqg_short_course"
library(predictionet)


###################################################
### chunk number 11: foopredictionet eval=FALSE
###################################################
## #line 335 "predictionet_pqg_short_course"
## library(help=predictionet)


###################################################
### chunk number 12: helpnetinf eval=FALSE
###################################################
## #line 361 "predictionet_pqg_short_course"
## help(netinf)


###################################################
### chunk number 13: firstnetinf
###################################################
#line 367 "predictionet_pqg_short_course"
mynet <- netinf(data=data1, priors=priors.count, priors.count=TRUE, priors.weight=0.5, maxparents=3, method="regrnet", seed=54321)


###################################################
### chunk number 14: regrnetopofig eval=FALSE
###################################################
## #line 376 "predictionet_pqg_short_course"
## ## network topology
## mytopoglobal <- net2topo(net=mynet, coefficients=FALSE)
## library(network)
## xnet <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
## mynetlayout <- plot.network(x=xnet, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5)


###################################################
### chunk number 15: regrnetopofig
###################################################
#line 386 "predictionet_pqg_short_course"
## network topology
mytopoglobal <- net2topo(net=mynet, coefficients=FALSE)
library(network)
xnet <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
mynetlayout <- plot.network(x=xnet, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.6)


###################################################
### chunk number 16: edgecoldiff
###################################################
#line 404 "predictionet_pqg_short_course"
## preparing colors
mycols <- c("blue", "green", "red")
names(mycols) <- c("prior", "data", "both")
myedgecols <- matrix("white", nrow=nrow(mytopoglobal), ncol=ncol(mytopoglobal), dimnames=dimnames(mytopoglobal))
## prior only
myedgecols[mytopoglobal == 0 & priors.count >= 1] <- mycols["prior"]
## data only
myedgecols[mytopoglobal == 1 & priors.count < 1] <- mycols["data"]
## both in priors and data
myedgecols[mytopoglobal == 1 & priors.count >= 1] <- mycols["both"]


###################################################
### chunk number 17: edgecoldiffig eval=FALSE
###################################################
## #line 417 "predictionet_pqg_short_course"
## mytopoglobal2 <- (myedgecols != "white") + 0
## ## network topology
## xnet2 <- network(x=mytopoglobal2, matrix.type="adjacency", directed=TRUE, loops=TRUE, vertex.attrnames=dimnames(mytopoglobal2)[[1]])
## plot.network(x=xnet2, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col=myedgecols, coord=mynetlayout)


###################################################
### chunk number 18: edgecoldiffig
###################################################
#line 426 "predictionet_pqg_short_course"
mytopoglobal2 <- (myedgecols != "white") + 0
## network topology
xnet2 <- network(x=mytopoglobal2, matrix.type="adjacency", directed=TRUE, loops=TRUE, vertex.attrnames=dimnames(mytopoglobal2)[[1]])
plot.network(x=xnet2, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col=myedgecols, coord=mynetlayout)


###################################################
### chunk number 19: regrnetinfcv
###################################################
#line 447 "predictionet_pqg_short_course"
myres.cv <- netinf.cv(data=data1, categories=3, priors=priors.count, maxparents=3, priors.weight=0.5, method="regrnet", nfold=10, seed=54321)


###################################################
### chunk number 20: rescv
###################################################
#line 453 "predictionet_pqg_short_course"
print(str(myres.cv, 1))


###################################################
### chunk number 21: edgecol
###################################################
#line 478 "predictionet_pqg_short_course"
## preparing colors
ii <- 0:10
mycols <- c("#BEBEBE", rev(rainbow(10, v=0.8, alpha=0.5)))
names(mycols) <- as.character(ii/10)
myedgecols <- matrix("#00000000", nrow=nrow(mytopoglobal), ncol=ncol(mytopoglobal), dimnames=dimnames(mytopoglobal))
for(k in 1:length(mycols)) { myedgecols[myres.cv$edge.stability == names(mycols)[k]] <- mycols[k] }
myedgecols[!mytopoglobal] <- "#00000000"


###################################################
### chunk number 22: edgestab eval=FALSE
###################################################
## #line 488 "predictionet_pqg_short_course"
## def.par <- par(no.readonly=TRUE)
## layout(rbind(1,2), heights=rbind(8,1))
## par(mar=c(1,1,1,1))
## ## network topology
## xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
## plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col=myedgecols, coord=mynetlayout)
## par(mar=c(0,3,1,3))
## plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="Stability scale", cex.main=1)
## rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
## text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
## par(def.par)


###################################################
### chunk number 23: edgestabfig
###################################################
#line 504 "predictionet_pqg_short_course"
def.par <- par(no.readonly=TRUE)
layout(rbind(1,2), heights=rbind(8,1))
par(mar=c(1,1,1,1))
## network topology
xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col=myedgecols, coord=mynetlayout)
par(mar=c(0,3,1,3))
plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="Stability scale", cex.main=1)
rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
par(def.par)


###################################################
### chunk number 24: genecolr2
###################################################
#line 536 "predictionet_pqg_short_course"
myr2 <- apply(myres.cv$prediction.score.cv$r2, 2, mean, na.rm=TRUE)
myr2 <- round(myr2*10)/10
## preparing colors
ii <- 0:10
mycols <- c("#BEBEBE", rev(rainbow(10, v=0.8, alpha=0.5)))
names(mycols) <- as.character(ii/10)
myvertexcols <- mycols[match(myr2, names(mycols))]
names(myvertexcols) <- names(myr2)


###################################################
### chunk number 25: edgestab eval=FALSE
###################################################
## #line 547 "predictionet_pqg_short_course"
## def.par <- par(no.readonly=TRUE)
## layout(rbind(1,2), heights=rbind(8,1))
## par(mar=c(1,1,1,1))
## ## network topology
## xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
## plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col=myvertexcols, jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout)
## par(mar=c(0,3,1,3))
## plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="R^2", cex.main=1)
## rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
## text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
## par(def.par)


###################################################
### chunk number 26: genepredabr2fig
###################################################
#line 563 "predictionet_pqg_short_course"
def.par <- par(no.readonly=TRUE)
layout(rbind(1,2), heights=rbind(8,1))
par(mar=c(1,1,1,1))
## network topology
xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col=myvertexcols, jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout)
par(mar=c(0,3,1,3))
plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="$R^2$", cex.main=1)
rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
par(def.par)


###################################################
### chunk number 27: genecolmcc
###################################################
#line 583 "predictionet_pqg_short_course"
mymcc <- apply(myres.cv$prediction.score.cv$mcc, 2, mean, na.rm=TRUE)
mymcc <- round(mymcc*20)/20
mymcc[mymcc < 0.5] <- 0.5
## preparing colors
ii <- 0:10
mycols <- c("#BEBEBE", rev(rainbow(10, v=0.8, alpha=0.5)))
names(mycols) <- seq(0.5, 1, 0.05)
myvertexcols <- mycols[match(mymcc, names(mycols))]
names(myvertexcols) <- names(mymcc)


###################################################
### chunk number 28: genepredabmccfig
###################################################
#line 597 "predictionet_pqg_short_course"
def.par <- par(no.readonly=TRUE)
layout(rbind(1,2), heights=rbind(8,1))
par(mar=c(1,1,1,1))
## network topology
xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col=myvertexcols, jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout)
par(mar=c(0,3,1,3))
plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="$MCC$", cex.main=1)
rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
par(def.par)


###################################################
### chunk number 29: validationpred
###################################################
#line 625 "predictionet_pqg_short_course"
## compute predictions
mynet.test.pred <-  netinf.predict(net=mynet, data=data2)
## performance estimation: R2
mynet.test.r2 <- pred.score(data=data2, pred=mynet.test.pred, method="r2")
## performance estimation: MCC
mynet.test.mcc <- pred.score(data=data2, categories=3, pred=mynet.test.pred, method="mcc")


###################################################
### chunk number 30: genecolr2test
###################################################
#line 635 "predictionet_pqg_short_course"
mynet.test.r2 <- round(mynet.test.r2*10)/10
## preparing colors
ii <- 0:10
mycols <- c("#BEBEBE", rev(rainbow(10, v=0.8, alpha=0.5)))
names(mycols) <- as.character(ii/10)
myvertexcols <- mycols[match(mynet.test.r2, names(mycols))]
names(myvertexcols) <- names(mynet.test.r2)


###################################################
### chunk number 31: genepredabr2testfig
###################################################
#line 647 "predictionet_pqg_short_course"
def.par <- par(no.readonly=TRUE)
layout(rbind(1,2), heights=rbind(8,1))
par(mar=c(1,1,1,1))
## network topology
xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col=myvertexcols, jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout)
par(mar=c(0,3,1,3))
plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="$R^2$", cex.main=1)
rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
par(def.par)


###################################################
### chunk number 32: genecolmcctest
###################################################
#line 665 "predictionet_pqg_short_course"
mynet.test.mcc <- round(mynet.test.mcc*20)/20
mynet.test.mcc[mynet.test.mcc < 0.5] <- 0.5
## preparing colors
ii <- 0:10
mycols <- c("#BEBEBE", rev(rainbow(10, v=0.8, alpha=0.5)))
names(mycols) <- seq(0.5, 1, 0.05)
myvertexcols <- mycols[match(mynet.test.mcc, names(mycols))]
names(myvertexcols) <- names(mynet.test.mcc)


###################################################
### chunk number 33: genepredabmcctestfig
###################################################
#line 678 "predictionet_pqg_short_course"
def.par <- par(no.readonly=TRUE)
layout(rbind(1,2), heights=rbind(8,1))
par(mar=c(1,1,1,1))
## network topology
xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col=myvertexcols, jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout)
par(mar=c(0,3,1,3))
plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="$MCC$", cex.main=1)
rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
par(def.par)


###################################################
### chunk number 34: regnetcvpriorsweight
###################################################
#line 705 "predictionet_pqg_short_course"
## priors.weight 0
myres.cv.pw0 <- netinf.cv(data=data1, categories=3, priors=priors.count, maxparents=3, priors.weight=0, method="regrnet", nfold=10, seed=54321)
## priors.weight 0.25
myres.cv.pw025 <- netinf.cv(data=data1, categories=3, priors=priors.count, maxparents=3, priors.weight=0.25, method="regrnet", nfold=10, seed=54321)
## priors.weight 0.5
myres.cv.pw050 <- myres.cv
## priors.weight 0.75
myres.cv.pw075 <- netinf.cv(data=data1, categories=3, priors=priors.count, maxparents=3, priors.weight=0.75, method="regrnet", nfold=10, seed=54321)
## priors.weight 0
myres.cv.pw1 <- netinf.cv(data=data1, categories=3, priors=priors.count, maxparents=3, priors.weight=1, method="regrnet", nfold=10, seed=54321)


###################################################
### chunk number 35: regnetcvpriorsweightfig eval=FALSE
###################################################
## #line 721 "predictionet_pqg_short_course"
## def.par <- par(no.readonly=TRUE)
## layout(mat=matrix(1:4, nrow=2, ncol=2, byrow=TRUE))
## ## priors.weight 0
## mytopot <- net2topo(net=myres.cv.pw0$net, coefficients=FALSE)
## xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
## plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 0 (data only)")
## ## priors.weight 0.25
## mytopot <- net2topo(net=myres.cv.pw025$net, coefficients=FALSE)
## xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
## plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 0.25")
## ## priors.weight 0.75
## mytopot <- net2topo(net=myres.cv.pw075$net, coefficients=FALSE)
## xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
## plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 0.75")
## ## priors.weight 1
## mytopot <- net2topo(net=myres.cv.pw1$net, coefficients=FALSE)
## xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
## plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 1 (priors only)")
## par(def.par)


###################################################
### chunk number 36: regnetcvpriorsweightfig2
###################################################
#line 745 "predictionet_pqg_short_course"
def.par <- par(no.readonly=TRUE)
layout(mat=matrix(1:4, nrow=2, ncol=2, byrow=TRUE))
## priors.weight 0
mytopot <- net2topo(net=myres.cv.pw0$net, coefficients=FALSE)
xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 0 (data only)")
## priors.weight 0.25
mytopot <- net2topo(net=myres.cv.pw025$net, coefficients=FALSE)
xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 0.25")
## priors.weight 0.75
mytopot <- net2topo(net=myres.cv.pw075$net, coefficients=FALSE)
xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 0.75")
## priors.weight 1
mytopot <- net2topo(net=myres.cv.pw1$net, coefficients=FALSE)
xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 1 (priors only)")
par(def.par)


###################################################
### chunk number 37: boxplotstabpwcvfig eval=FALSE
###################################################
## #line 776 "predictionet_pqg_short_course"
## gg <- c("0", "0.25", "0.50", "0.75", "1")
## mystab.cv.pw <- list(myres.cv.pw0$edge.stability[myres.cv.pw0$topology == 1], myres.cv.pw025$edge.stability[myres.cv.pw025$topology == 1], myres.cv.pw050$edge.stability[myres.cv.pw050$topology == 1], myres.cv.pw075$edge.stability[myres.cv.pw075$topology == 1], myres.cv.pw1$edge.stability[myres.cv.pw1$topology == 1])
## names(mystab.cv.pw) <- gg
## boxplot(mystab.cv.pw, xlab="priors.weight", ylim=c(0, 1), ylab="Edge stability", border="grey", col="lightgrey")
## points(x=jitter(x=rep(1:length(mystab.cv.pw), times=sapply(mystab.cv.pw, length)), amount=0.15), y=unlist(mystab.cv.pw), cex=0.55, pch=16, col="royalblue")


###################################################
### chunk number 38: boxplotstabpwcvfig2
###################################################
#line 786 "predictionet_pqg_short_course"
gg <- c("0", "0.25", "0.50", "0.75", "1")
mystab.cv.pw <- list(myres.cv.pw0$edge.stability[myres.cv.pw0$topology == 1], myres.cv.pw025$edge.stability[myres.cv.pw025$topology == 1], myres.cv.pw050$edge.stability[myres.cv.pw050$topology == 1], myres.cv.pw075$edge.stability[myres.cv.pw075$topology == 1], myres.cv.pw1$edge.stability[myres.cv.pw1$topology == 1])
names(mystab.cv.pw) <- gg
boxplot(mystab.cv.pw, xlab="priors.weight", ylim=c(0, 1), ylab="Edge stability", border="grey", col="lightgrey", outline=FALSE)
points(x=jitter(x=rep(1:length(mystab.cv.pw), times=sapply(mystab.cv.pw, length)), amount=0.15), y=unlist(mystab.cv.pw), cex=0.55, pch=16, col="royalblue")


###################################################
### chunk number 39: boxplotr2pwcvfig eval=FALSE
###################################################
## #line 801 "predictionet_pqg_short_course"
## gg <- c("0", "0.25", "0.50", "0.75", "1")
## myr2.cv.pw <- cbind(apply(myres.cv.pw0$prediction.score.cv$r2, 2, mean, na.rm=TRUE), apply(myres.cv.pw025$prediction.score.cv$r2, 2, mean, na.rm=TRUE), apply(myres.cv.pw050$prediction.score.cv$r2, 2, mean, na.rm=TRUE), apply(myres.cv.pw075$prediction.score.cv$r2, 2, mean, na.rm=TRUE), apply(myres.cv.pw1$prediction.score.cv$r2, 2, mean, na.rm=TRUE))
## colnames(myr2.cv.pw) <- gg
## gg <- factor(rep(gg, each=genen), levels=gg)
## boxplot(myr2.cv.pw, xlab="priors.weight", ylim=c(0, 1), ylab="$R^2$", border="grey", col="lightgrey", outline=FALSE)
## points(x=jitter(x=rep(1:ncol(myr2.cv.pw), times=nrow(myr2.cv.pw)), amount=0.15), y=as.vector(t(myr2.cv.pw)), cex=0.55, pch=16, col="royalblue")


###################################################
### chunk number 40: boxplotr2pwcvfig2
###################################################
#line 812 "predictionet_pqg_short_course"
gg <- c("0", "0.25", "0.50", "0.75", "1")
myr2.cv.pw <- cbind(apply(myres.cv.pw0$prediction.score.cv$r2, 2, mean, na.rm=TRUE), apply(myres.cv.pw025$prediction.score.cv$r2, 2, mean, na.rm=TRUE), apply(myres.cv.pw050$prediction.score.cv$r2, 2, mean, na.rm=TRUE), apply(myres.cv.pw075$prediction.score.cv$r2, 2, mean, na.rm=TRUE), apply(myres.cv.pw1$prediction.score.cv$r2, 2, mean, na.rm=TRUE))
colnames(myr2.cv.pw) <- gg
gg <- factor(rep(gg, each=genen), levels=gg)
boxplot(myr2.cv.pw, xlab="priors.weight", ylim=c(0, 1), ylab="$R^2$", border="grey", col="lightgrey", outline=FALSE)
points(x=jitter(x=rep(1:ncol(myr2.cv.pw), times=nrow(myr2.cv.pw)), amount=0.15), y=as.vector(t(myr2.cv.pw)), cex=0.55, pch=16, col="royalblue")


###################################################
### chunk number 41: regnetcvpriorsweightcompr2
###################################################
#line 826 "predictionet_pqg_short_course"
## Friedman test to test whether at least one of the networks gives statically different predictive ability
print(friedman.test(y=myr2.cv.pw))

## Pairwise Wilcoxon Rank Sum tests
print(pairwise.wilcox.test(x=as.vector(myr2.cv.pw), g=gg, paired=TRUE, exact=FALSE, alternative="two.sided", p.adjust.method="bonferroni"))


###################################################
### chunk number 42: boxplotr2pwtestfig2
###################################################
#line 839 "predictionet_pqg_short_course"
gg <- c("0", "0.25", "0.50", "0.75", "1")
myr2.test.pw <- cbind(pred.score(data=data2, pred=netinf.predict(net=myres.cv.pw0$net, data=data2), method="r2"), pred.score(data=data2, pred=netinf.predict(net=myres.cv.pw025$net, data=data2), method="r2"), pred.score(data=data2, pred=netinf.predict(net=myres.cv.pw050$net, data=data2), method="r2"), pred.score(data=data2, pred=netinf.predict(net=myres.cv.pw075$net, data=data2), method="r2"), pred.score(data=data2, pred=netinf.predict(net=myres.cv.pw1$net, data=data2), method="r2"))
colnames(myr2.test.pw) <- gg
gg <- factor(rep(gg, each=genen), levels=gg)
boxplot(myr2.test.pw, xlab="priors.weight", ylim=c(0, 1), ylab="$R^2$", border="grey", col="lightgrey", outline=FALSE)
points(x=jitter(x=rep(1:ncol(myr2.test.pw), times=nrow(myr2.test.pw)), amount=0.15), y=as.vector(t(myr2.test.pw)), cex=0.55, pch=16, col="royalblue")


###################################################
### chunk number 43: boxplotr2pwtestfig2
###################################################
#line 850 "predictionet_pqg_short_course"
gg <- c("0", "0.25", "0.50", "0.75", "1")
myr2.test.pw <- cbind(pred.score(data=data2, pred=netinf.predict(net=myres.cv.pw0$net, data=data2), method="r2"), pred.score(data=data2, pred=netinf.predict(net=myres.cv.pw025$net, data=data2), method="r2"), pred.score(data=data2, pred=netinf.predict(net=myres.cv.pw050$net, data=data2), method="r2"), pred.score(data=data2, pred=netinf.predict(net=myres.cv.pw075$net, data=data2), method="r2"), pred.score(data=data2, pred=netinf.predict(net=myres.cv.pw1$net, data=data2), method="r2"))
colnames(myr2.test.pw) <- gg
gg <- factor(rep(gg, each=genen), levels=gg)
boxplot(myr2.test.pw, xlab="priors.weight", ylim=c(0, 1), ylab="$R^2$", border="grey", col="lightgrey", outline=FALSE)
points(x=jitter(x=rep(1:ncol(myr2.test.pw), times=nrow(myr2.test.pw)), amount=0.15), y=as.vector(t(myr2.test.pw)), cex=0.55, pch=16, col="royalblue")


###################################################
### chunk number 44: regnettestpriorsweightcompr2
###################################################
#line 862 "predictionet_pqg_short_course"
## Friedman test to test whether at least one of the networks gives statically different predictive ability
print(friedman.test(y=myr2.test.pw))

## Pairwise Wilcoxon Rank Sum tests
print(pairwise.wilcox.test(x=as.vector(myr2.test.pw), g=gg, paired=TRUE, exact=FALSE, alternative="two.sided", p.adjust.method="bonferroni"))


###################################################
### chunk number 45: sessioninfo
###################################################
#line 892 "predictionet_pqg_short_course"
toLatex(sessionInfo(), locale=FALSE)


